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Great attention is given to the first star formation and the epoch of reionization as main targets 
of planned large radio interferometries (e.g. Square Kilometre Array). Recently, it is claimed that 
the supersonic relative velocity between baryons and cold dark matter can suppress the abundance 
of first stars and impact the cosmological reionization process. Therefore, in order to compare 
observed results with theoretical predictions it is important to examine the effect of the supersonic 
relative motion on the small-scale structure formation. In this paper, we investigate this effect on the 
nonlinear structure formation in the context of the spherical collapse model in order to understand 
the fundamental physics in a simple configuration. We show the evolution of the dark matter sphere 
with the relative velocity by both using N-body simulations and numerically calculating the equation 
of motion for the dark matter mass shell. The effects of the relative motion in the spherical collapse 
model appear as the delay of the collapse time of dark matter halos and the decrease of the baryon 
mass fraction within the dark matter sphere. Based on these results, we provide the fitting formula 
of the critical density contrast for collapses with the relative motion effect and calculate the mass 
function of dark matter halos in the Press-Schechter formalism. As a result, the relative velocity 
decreases the abundance of dark matter halos whose mass is smaller than 10^ Mq/Zi. 


I. INTRODUCTION 

The standard cosmological model, called the ACDM 
model, composed with two relativistic species (photons 
and neutrinos), two nonrelativistic matters (baryons and 
dark matter), and the energy having negative pres¬ 
sure (dark energy) with a nearly scale-invariant spec¬ 
trum of curvature perturbations, has achieved great suc¬ 
cess in explaining large-scale cosmological observations, 
e.g., large-scale structure formation [i| and cosmic mi¬ 
crowave background [2[ . 

The theoretical description of small-scale structure for¬ 
mation is, however, still debatable. Understanding of 
small-scale structure formation at high redshifts is es¬ 
sential to study first stars and the epoch of reioniza¬ 
tion (EoR). One expects redshifted 21 cm lines from 
the hyperfine structure of hydrogen atoms as the pow¬ 
erful probe for the EoR and first stars d, Q and the 
matter density underlying HI distribution constrains ex¬ 
tended parameters of the ACDM model Currently, 

to probe such small-scale structure formation, there are 
many planed observations including Murchison Widefield 
Array [H and Square Kilometre Array (SKA) [9|. There¬ 
fore, nowadays, the detailed studies on small-scale struc¬ 
ture formation at high redshifts attract a lot of attention. 

Recently, Ref. 0 reports the importance of the super¬ 
sonic relative motion between dark matter and baryons 
on small-scale structure formation related to the EoR. 
This supersonic relative motion is originated from the 
difference in the motions between baryons and dark mat¬ 
ter before recombination. Baryons before recombination 
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are tightly coupled with photons by Thomson scattering. 
As a result, baryons and photons act as one fluid with 
the sound speed ^ c/\/3 and have the velocity field asso¬ 
ciated with the acoustic oscillation. On the other hand, 
dark matter does not suffer from Thomson scattering and 
dark matter density fluctuations can grow gravitation¬ 
ally. Therefore, the relative motion between baryons and 
dark matter is induced. After recombination, baryons 
are fully decoupled with photons and the sound speed of 
baryons quickly drops to ^ 6 km/s. Since the root mean 
square of the relative velocity reaches ^ 30 km/s at that 
time, the relative velocity is about five times larger than 
the sound speed of baryons. Because the relative motion 
is highly supersonic, the effect on the structure forma¬ 
tion could be significant. In particular, the abundance 
of small dark matter halos (M ^ 10^ ^ 0 ) is highly sup¬ 
pressed due to the supersonic relative motion. This ef¬ 
fect has been intensively studied by many authors with 
N-body/smoothed particle hydrodynamics (SPH) simu¬ 
lations [lll-[l^. Therefore, according to the effect on the 
structure formation, the scenario of the cosmic reioniza¬ 
tion and the prediction of the 21 cm line signals from the 
EoR could be modified from those predicted in the con¬ 
ventional cosmological model. The recent relevant stud¬ 
ies are reviewed in Ref. 0 . 

So far the effects of the relative motion on the structure 
formation have been studied in numerical simulations 
mainly, because these effects are complicated. However, 
a study with the analytical model is useful to obtain some 
insights into physics involved in complicated phenomena. 
Additionally the analysis of observation data with nu¬ 
merical simulations generally takes enormous time and, 
sometimes, it seems unrealistic. Therefore, modeling in a 
form which is easy to handle in analytic studies is highly 
required. 
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In the study on the structure formation, the halo mass 
function is one of the interesting quantities. In particular, 
the Press-Schechter formalism with the sphere collapse 
model provides the mass function in the analytical form 
which relatively agrees well with the results of N-body 
simulations. 

In this paper, we revisit the effect of the supersonic 
relative motion on small-scale structure formation in the 
context of the spherical collapse model by both N-body 
simulations and a semianalytical way. The effect of the 
relative motion on the spherical collapse can be repre¬ 
sented as the modification of the critical density con¬ 
trast. We propose a fitting formula of the critical density 
contrast, as a function of the amplitude of the relative 
motion, the halo mass and the initial density fluctuation 
within dark matter halos. We also apply this fitting for¬ 
mula to evaluate the mass function of small dark matter 
halos based on the Press-Schechter formalism. 

This paper is organized as follows. In Sec. [IT] we re¬ 
view the effect of supersonic relative motion on the per¬ 
turbation theory by taking into account the background 
velocity of baryons and construct the spherical collapse 
model with two components. In Sec. IIIII we describe the 
setup of our N-body simulation, and we show the results 
of the N-body simulations and check the reproducibility 
of the spherical collapse in Sec. IIVI Moreover we present 
the change of the collapse time by supersonic relative 
motion and the validity of the semianalytical model in¬ 
troduced in Sec.[IIl Section[Vlis devoted to the discussion 
of the relative motion effect on dark matter halos. We 
discuss the modification of the baryon fraction in a dark 
matter halo by the relative motion, and, providing a fit¬ 
ting formula of the modified collapse time (i.e. the critical 
density contrast for the collapse). We show the suppres¬ 
sion of the dark matter halo abundance around the EoR 
as an application of our results. Finally we summarize 
this paper in Sec. 


II. ANALYTICAL FORMALISM 

In this section, we show the effect of the supersonic 
relative motion between baryons and dark matter on the 
linear perturbation theory, and evaluate this effect on 
the nonlinear growth by adopting the spherical collapse 
model. 


A. Perturbation theory 

Since the supersonic relative motion between baryons 
and dark matter has been analytically studied in the 
moving-background perturbation theory (MBPT) la, 
we first make a brief review of the MBPT. The MBPT 
introduces the background peculiar velocity which corre¬ 
sponds to the relative velocity between baryons and dark 
matter, i.e. = vi^^. According to the energy mo¬ 

mentum conservation equation with homogeneous back¬ 


ground densities of baryons and dark matter, the evolu¬ 
tion of Vbc is in reverse proportion to the scale factor due 
to the cosmic expansion. 

In the MBPT, the first order energy momentum con¬ 
servation equations after recombination are given by 
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where Cs is the sound velocity of the baryon fluid, the sub¬ 
scripts c and b denote cold dark matter and baryons re¬ 
spectively and 0 = ia~^V • v represents the divergence of 
the peculiar velocity. In Eq. m we take the frame where 
the background velocity of cold dark matter is absent. 
In other words, = Vhc and = 0. For simplicity, 
we ignore perturbations of the sound velocity although 
they might affect the growth of the density fluctuation 
on small scales [iM3- Equation ([T]) can be rewritten 
to the second order differential equations of the density 
fluctuations as 


^rrdSc^3H\ 

dt^ ~ ^ dt^ 2 + 

-^=-(2H + + —{^cdc + 

- (Cs + (2) 

where /i = • k/\vhc\\k\. 

Equation ^ tells us that the relative motion prevents 
the growth of density fluctuations on small scales in the 
same way of the fluid pressure in the discussion of the 
Jeans instability. On large scales where the relative mo¬ 
tion does not have the preferred direction, while the odd 
term of /i in the last equation of Eq. (|2]) vanishes by av¬ 
eraging over all random directions of the relative motion, 
the third term of the right-hand side is enhanced due to 
the existence of the term with As a result, the ef¬ 
fective Jeans scale (the suppression scale) of Eq. (|2]) be¬ 
comes large due to the existence of the relative velocity. 
Since the relative velocity after recombination is roughly 
r\j 5cs, the suppression scale for the relative mo¬ 
tion is khc = aUj^ 40 hMpc“^. The correspond¬ 
ing mass scale for the suppression is Mhc ^ 10^ M^jh. 

However, when we consider a local sufficiently small 
patch, the odd term of /i cannot vanish in the patch. 
Instead, the relative motion in this patch can be assumed 
to be a homogeneous flow with one direction. In this case, 
when the relative velocity is larger than the Hubble flow, 
liVhc > k/aH^ the density fluctuations inside the patch 
start to grow exponentially due to the relative motion 
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flow as shown in Eq. ©• The perturbation theory is not 
valid in this case, and we need to consider the effect of 
the relative velocity on the nonlinear growth, e.g., in the 
spherical collapse model. 


B. Spherical collapse model 

The spherical collapse model is a simple analytical 
model to investigate the nonlinear evolution of an over¬ 
density region. In this model, the evolution of the over¬ 
density region is described as the motion of the constant 
density spheres. Let us consider the collapse of a mass 
shell inside which the mass is M = 47rxfpi(l -1- ^i)/3, 
where Xi is the initial radius and 6i is the initial density 
contrast within the sphere with the radius Xi (hereafter 
the subscript i represents the initial time value). The 
equation of motion (EoM) for the proper radius x of a 
shell is written as 


d‘^x _ GM 

Equation (|3|) can be solved analytically, and the solution 
is given by 
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where ti is the initial time and Vi is the initial velocity. 
The solution of Eq. (j4]) depends only on 6i and vi. In 
order to keep the constant mass dM/dt = 0, we give the 
initial velocity of the shells 


Vi = HiXi 


3(1 + 5i) 


( 5 ) 


where we assume the matter dominated era, t oc 
The first term represents the Hubble flow and the second 
term corresponds to the peculiar velocity. According to 
Eq. ®, the solution, Eq. ® depends on only the initial 
density fluctuation, 5i. Eurthermore, we can obtain the 
critical density contrast that is the density contrast at the 
collapse time ^ = 27r in the linear perturbation theory. 


^ _a{e = 2 tt) ^ 
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where we use the fact that the growth factor of the matter 
density perturbation is proportional to the scale factor in 
the matter dominated era. 

Next we consider the effect of the homogeneous su¬ 
personic relative motion on the nonlinear evolution in 
the spherical collapse model. A simple extension is to 
introduce the two kinds of mass shells for dark matter 
and baryons. Taking into account the supersonic relative 


motion, the baryon mass shells have the initial bulk ve¬ 
locity, because we take the frame where baryons have the 
homogeneous relative flow to dark matter. As a result, 
the collapsing of the baryon mass shells is not spheri¬ 
cal and these mass shells are collapsing to the different 
position from the dark matter shells. However, the col¬ 
lapse of dark matter precedes the one of the baryons in a 
halo formation and we are interested in a baryon fraction 
within a collapsed dark matter halo. Therefore, we focus 
on the dark matter mass shell and, instead of following 
the evolution of the baryon mass shells, we introduce the 
baryon mass within the dark matter mass shell, M^, and 
rewrite Eq. (j3]) to 

d^xc _ G{M,^i + M,) 

dt‘^ x‘1 

Mc,i — Pc,i^c,ii^ ^c,i)’! 

Anr 

Mh=—pbxlil + db), (7) 

where Mc^i is the mass within the shell with the initial 
radius at Xc is the radius of the mass shell with Mc^i 
at each time, and is the baryon mass in the dark 
matter shell at the radius Xc. As a shell collapses, the 
baryon mass inside the shell increases with the growth 
of 6b, namely the baryon collapsing. Although the den¬ 
sity fluctuation of baryons without the relative motion 
catches up soon with that of dark matter, the relative 
motion prevents this process. Therefore, the effect of the 
relative motion is included through the evolution of M^. 
However it is difficult to evaluate analytically Mb with 
the relative motion. Thus in order to compute Eq. 0, 
we adopt Mb obtained from the N-body simulation in 
the following section. We also compare the result based 
on the spherical collapse model with that from the full 
N-body simulations in the later section. 


III. N-BODY SIMULATION 

Besides the analytical way mentioned in the previous 
section, we evaluate the effect of the supersonic motion 
between dark matter and baryons on the structure for¬ 
mation at high redshifts by using N-body simulations. 
In this section, we describe the setup of our N-body 
simulations. We perform N-body simulations with the 
public code Gadget-2 [T^. In all N-body simulations, 
the cosmological parameters are set to (O^, Oa h) = 
(0.31, 0.69, 0.68) with Vtb/^m ^ 1/6. The effect of 
the supersonic relative motion on the structure forma¬ 
tion works after the decoupling between photons and 
baryons. Therefore, the initial redshift for the simula¬ 
tions is Zi = 1000. Note that our results almost do not 
depend on the cosmological parameters because we are 
interested in the structure formation in the matter dom¬ 
inated era. 

Eor the initial distribution of the particles, we consider 
the spherical top-hat overdensity region in the isolated 
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system. The simulation box has the uniform distribution 
of the particles with a uniform overdensity sphere. We set 
the box size to Tbox = 200 kpc/h and the radius of the 
overdensity sphere to = 50 kpc/h. Note that we denote 
hereafter x as the proper distance and r as the comoving 
distance. In the box, the number of the uniform particles 
is 3 X 10^ and the initial density contrast of dark matter 
in the overdense sphere is Sc,i = 0.033. Thus the mass of 
particles is 2 x 10^ Mq/H and the mass of dark matter 
within the initial overdensity sphere is given by Me ^ 
4 X 10^ MqIH. Moreover we set the softening parameter 
to e = 0.1 kpc/h. We confirm that changes in these 
parameters do not affect our result qualitatively. 

According to the cosmological perturbation theory, the 
amplitude of baryon density fluctuations is 1% of that of 
dark matter density fluctuations with k = 100 Mpc“^ at 
z ^ 1000. The initial fluctuations of baryons are negli¬ 
gible compared with those of dark matter at Zi = 1000. 
Therefore, we assume that 1/6) of the uniform 

distributed particles is composed of baryons and there is 
no baryon fluctuation in the over density sphere. 

Figure [T] shows the initial configuration of the particles. 
The red dots represent the particle uniformly distributed 
in the box and the green dots are for the particles in¬ 
cluded in the overdense sphere. Figure [2] shows the initial 
density contrast of dark matter particles as a function of 
the radius from the center. In this figure, the red points 
indicate the values averaged over five realizations of our 
simulations and the error bars represent the shot noise 
caused by the finite particle number. The black line is 
the analytical prediction from our initial condition. This 
figure tells us that the density is constant in the top-hat 
sphere. Therefore we can convert the initial position of 
mass shell to the mass contained within each shell by us¬ 
ing the relation Mc^i = 47rp(l+(i^)x?/3. We set the initial 
velocity of dark matter given by only the second term of 
Eq. dSj), because N-body simulations are performed in the 
comoving coordinate. 

In order to take into account the supersonic relative 
motion, we give the additional velocity to all baryons. 
The correlation of the supersonic relative velocity has the 
significant value on larger scale than scales of our interest 
that are smaller than Mpc. Therefore, we assume that all 
baryons in the simulations have the constant supersonic 
relative velocity in one direction. In other words, in 
the simulation, the additional initial velocity for baryons 
is represented as = (u^c^O^O)- 

All terms related to the relative velocity in Eq. © are 
proportional to kv^jc- This fact suggests that the effect of 
relative motion on the spherical collapse of dark matter 
halos also depends on a factor v^ck oc vijc/mI^^ . Thus, in¬ 
stead of changing both the dark matter halo mass Me and 
the relative velocity Ubc, we perform numerical simula¬ 
tions for different relative velocities (5 km, 15 km, 30 km, 
50 km, 100 km, 150 km, 200 km, 300 km, 500 km) with 
fixing the dark matter halo mass Me ^ 4 x 10^ Mq/Zi, 
in order to evaluate the dependence of the effect of the 
relative velocity on Me and v^e- 
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FIG. 1: The initial configuration of the particles. The green 
particles are contained within the top-hat sphere and the red 
points are otherwise. 



FIG. 2: The initial density contrast distribution of dark mat¬ 
ter. The red points are the values averaged five realizations 
and the error bars show the shot noise caused by the number 
of particles contained within mass shells. The black line is 
the analytical prediction. 


In the simulations, we use the periodic boundary con¬ 
dition. Therefore, when the relative velocity is higher 
than 200 km/s, baryon particles leaving the simulation 
box along the direction of the relative velocity reenter the 
box from the opposite direction due to the boundary con¬ 
dition. In this case, the distribution of reentering baryons 
is no longer homogeneous in the perpendicular direction 
to the relative velocity and, resultantly, the dependence 
on the boundary condition arises in the results. In order 
to remove this dependence, we make the perpendicular 
positions of baryons random when the baryon particles 
reenter the simulation box. 
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FIG. 3: The evolutions of radii of mass shells containing 
0.6Mc (red), O.SMc (green), and Me (blue). The solid black 
line is the analytical solution of spherical collapse model with 
— 0.028. The dashed black line is the analytical solution 
corrected in consideration of the periodic boundary condition 
with B — A. 


FIG. 4: The turnaround time (upper panel) and the reference 
collapse time (lower panel). The red points are the results 
of N-body simulations with the standard error measured five 
realizations, and the black lines are the analytical predictions. 
In the upper panel the shaded region shows the error caused 
from shot noise shown in Fig. The dashed line shows the 
prediction from the solutions of Eq. dS]) with B converted from 

Ti. 


IV. RESULT 

In this section, we present the results of our N-body 
simulations. First we show the result of the reference 
model that is the case without the relative velocity. 
Figure [3] shows the time evolutions of the radii of mass 
shells from the center. In our simulation, we determined 
the center of the collapsing shells by using the mean posi¬ 
tion of the particles contained initially within the top-hat 
over density sphere at each time step. In this figure, the 
red line represents the mass shell containing 0.6Mc (r^ = 
42 kpc/h), and the green and blue lines are for that con¬ 
taining O.SMc {vi = 46 kpc/h) and Me (n = 50 kpc/h), 
respectively. Additionally the black line corresponds to 
the analytical solution of the spherical collapse model, 
Eq. dH), with Sm,i = (Pc,i^c,i+P6,i^6,i)/(Pc,i+P6,i) = 0.028 
and the velocity Vi given by Eq. (jHj). Note that the shell 
evolution of the spherical collapse model depends on 6m,i 
only. Therefore, in our initial condition where the over¬ 
density sphere has the homogeneous density profile, the 
spherical collapse model predicts that all mass shells in¬ 
side the overdensity sphere trace the black dashed line 
and collapse at the same time, independently on the mass 
contained by the mass shells. 

The evolutions of radii of the mass shells from N-body 
simulations agree with the analytic solution of the spher¬ 
ical collapse model before the turnaround time when 
the radius reaches the maximum. However, after the 
turnaround time, the results from N-body simulations 
deviate from the analytic solution. One of the reasons 
for this deviation is that the particles cannot be concen¬ 
trated on the infinitesimal point in N-body simulations. 
Therefore, the particles in the simulations begin to be 
relaxed with each other after the turnaround time and 


the collapse is prevented. Eurt her more the shot noise in¬ 
duces the substructures inside the collapsing sphere and 
the ejection of the particles from the mass shells, and 
resultantly causes the dispersion of the collapse time as 
discussed in Ref. [T^ . These effects of relaxation and shot 
noise lead the delay of the collapse and cause the devi¬ 
ation from the analytical solution after the turnaround 
time. 

Eigure[3]also shows that the outer mass shells collapse 
later than the inner ones, although the theoretical spher¬ 
ical collapse model claims that all mass shells collapse 
at the same time. The reason is the effect of the peri¬ 
odic boundary condition. In order to evaluate this effect 
simply, we consider the motion of a particle along the 
x-axis direction with the periodic boundary condition. 
Since we should take into account the gravitational force 
from the overdensity region in the other boxes due to the 
boundary condition, the EoM, Eq. m, along the x-axis 
is corrected to 


(Px 

dP 


2(1 + Si) 2Si 

_9(naB — x)^ 


2Si 

9(naB + x)^ _ 
( 8 ) 


where B = aiLBo^jxi. When the position of the shell 
is close to the center of the simulation box at the ini¬ 
tial time, B becomes small and vice versa. Namely, the 
smaller B is, the more efficient the boundary effect is. 
In our calculation, the most outer shell that contains the 
mass Me inside has B = A. The evolution for H = 4 is 
plotted in the black dashed line in Eig. [3l 

EigureS] shows the turnaround time (upper panel) and 
the reference collapse time (lower panel) in the reference 
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FIG. 5: The evolutions of mass shell containing 0.6Mc with 
supersonic relative velocities Vbc = 30 km (red), Vhc = 50 km 
(green), Vbc = 100 km (blue), Vbc = 150 km/s (magenta) and 
reference model (black). Each dashed line is the solution of 
Eq. © with the baryon density fluctuation derived from N- 
body simulations. 


model as functions of the initial radius of the mass shell. 
Here, the turnaround time and the reference collapse time 
are defined as the times when the radius of the mass shell 
becomes maximum and minimum, respectively. In this 
figure, the red points with the error bars represent the 
averages and the standard errors obtained from the five 
realizations of N-body simulations, and the black line 
corresponds to the theoretical predictions in the spheri¬ 
cal collapse model. Moreover the dashed line shows the 
turnaround time obtained from Eq. m which includes 
the effect of the boundary condition. One can find that 
the turnaround time is consistent with the theoretical 
prediction within = 40 kpc/h. The one of the reasons 
why the outer mass shells turn around later is the ef¬ 
fect of the periodic boundary condition discussed above. 
In the case where we take the boundary condition into 
account, the turnaround time matches the theoretical 
prediction within = 43 kpc/h that corresponds with 
M = 0.65Mc. The difference between the turnaround 
time estimated from the outer shells than Vi = 43 kpc/h 
and theoretical one is due to the ejections of particles. As 
we have mentioned, the reference collapse time in N-body 
simulations delays for all , compared with the theoreti¬ 
cal prediction. Additionally, similarly to the turnaround 
time, the deviation becomes large as Vi increases. 

Next we show how the supersonic relative motion af¬ 
fects the collapse in N-body simulations. Performing five 
realizations for different we obtain the typical evo¬ 
lution of mass shells by averaging each realization. Fig¬ 
ure [5] represents the evolutions of the mass shells contain¬ 
ing 0.6Mc with four different supersonic relative veloci¬ 
ties, Vbc = 30 km/s (in red), Vbc = 50 km/s (in green), 
Vbc = 100 km/s (in blue) and Vbc = 150 km/s (in ma¬ 
genta). We can convert the results for different velocities 



FIG. 6: The turnaround time (upper panel) and the reference 
collapse time (lower panel) as the function of the initial radius 
of the mass shell with the supersonic relative velocity Vbc = 
30 km (red), Vbc = 50 km (green), Vbc = 100 km (blue), 
Vbc = 150 km/s (magenta) and the reference model (black). 
In the upper panel the dashed lines are the turnaround time 
estimated from the solution of the semianalytical model. 


with a fixed mass into those for different masses with 
a fixed velocity through the dependence of the relative 
velocity effect on Vbd^c as mentioned above. 

For comparison, the corresponding evolution of the 
mass shell with 0.6Mc in the reference model {vbc — 
0 km/s) is plotted as the black solid line. As the rel¬ 
ative velocity becomes large, the start of the collapse 
delays and the maximum radius increases. Therefore, 
we can conclude that the supersonic relative motion pre¬ 
vents the collapse. Additionally we show the solutions 
of Eq. ([7j) as the dashed lines in Fig. [5l Solving Eq. ([7]) 
numerically, we use the baryon fluctuation d obtained 
from the particle data of the N-body simulations within 
Ti < 42 kpc/h corresponding to 0.6Mc. We find that the 
semianalytical model agrees with the results of N-body 
simulations before the turnaround time. 

To illustrate the delay of the collapse due to the su¬ 
personic relative motion we plot the turnaround time 
and the reference collapse time in Fig. [6] for the differ¬ 
ent supersonic relative velocities. In this figure, both 
the turnaround and the reference collapse time are rep¬ 
resented as the functions of the initial radius of the mass 
shell. We show additionally the turnaround times ob¬ 
tained from the solutions of Eq. ([7|) as the dashed lines in 
the top panel of Fig.[6l The evaluations of the turnaround 
time from the semianalytical solutions are consistent with 
the turnaround times from N-body simulations within 
Ti ^ 40 kpc//i that is same as the reference case. In the 
following discussions, we use twice the turnaround time 
as the collapse time. In this case, the scale factor at the 

collapse time is given by Ucoi = (2tta)^^^ 
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metric distribution for baryons, the dark matter collapse 
remains spherical. Around the reference collapse time 
of the dark matter shells, the baryon fraction estimated 
within Ti = 42 kpc/h starts to oscillate. This is mainly 
due to the difference of the density peak positions be¬ 
tween baryons and dark matter. The baryon over density 
region is attracted by that of dark matter gravitationally 
and oscillates around. As time goes, the difference of 
the peak positions will be relaxed and the peak position 
of baryons is expected to overlap that of dark matter. 
Note that this result is based on the spherical collapse 
model which is an ideal isolated system. However, the 
actual collapse happens with many surrounding effects 
as shown in cosmological simulations. Therefore, to eval¬ 
uate the baryon fraction properly, these effects could be 
not negligible. 


FIG. 7: The baryon fractions within the dark matter over¬ 
density sphere whose initial radius are 50 kpc/h (solid lines) 
and 42 kpc/h (dashed lines). The arrows show the reference 
collapse times. 


V. DISCUSSION 

In this section, we show the baryon fraction within the 
dark matter over density sphere. We also discuss the de¬ 
lay of the collapse time and provide the fitting formula of 
the critical density contrast with the relative motion be¬ 
tween dark matter and baryons. Furthermore we present 
the modification of the halo mass function by taking ac¬ 
count of the relative motion. 


A. Baryon fraction 


B. Delay of the halo formation 

The supersonic relative motion delays the collapse time 
of dark matter halos as shown in Fig. [6l In the spherical 
collapse model without the relative velocity, the collapse 
time is dependent on only the initial density fluctuation. 
However, in the case of the nonzero relative velocity, the 
collapse time depends on the halo mass Me and the am¬ 
plitude of the relative velocity. Additionally, the effect 
of the relative motion cumulatively becomes large for the 
small initial density contrast, because the small initial 
density contrast takes longer time to the collapse. There¬ 
fore, the change of the collapse time with the relative ve¬ 
locity is represented as a function of Me, v^c and Sc,i We 
define the correction of the scale factor at the collapse 
time related with the modification of the critical density 
contrast, as 


First we consider the baryon fraction which repre¬ 
sents the mass ratio between baryons and total mat¬ 
ter (/^ = M\)/Mm) within the dark matter over-density 
sphere. The baryon fraction is important not only to es¬ 
timate the effect of the supersonic relative motion on the 
collapse of dark matter spheres, but also to discuss the 
first star formation or observables related with baryons. 
We calculate the baryon fraction by counting baryon 
particles within the dark matter collapsing over-density 
sphere. 

Figure [7| shows the time evolutions of the baryon frac¬ 
tion. As the relative velocity increases, the baryon frac¬ 
tion becomes smaller. In the case with the nonzero rela¬ 
tive velocity, the baryon overdensity region is no longer 
spherically symmetric and the peak position of baryon 
density is different from that of dark matter, depending 
on the amplitude of the relative velocity. However, such 
asymmetry of the baryon distribution does not affect the 
dark matter collapse well. As shown in Fig. [5l the dark 
matter collapse in the N-body simulations is consistent 
with the spherical collapse model until the turnaround 
time. Therefore, we can infer that, in spite of the asym- 


^c,i) 


^^col(Afc5^c,i) ^o(^c,i) 


(9) 


where ao((^c,i) is the scale factor at the collapse time with¬ 
out the relative velocity between baryons and dark mat¬ 
ter. 

Figure [8] shows the relative difference ^ as a func¬ 
tion of the relative velocity with a fixed mass Me ^ 
4x10^ Mq/ h. We plot A for different three initial density 
contrasts estimated from the N-body simulations with 
our boundary condition (shaded regions) or the usual pe¬ 
riodic boundary condition (dashed lines). We find that 
the effect of the supersonic relative velocity is negligible 
for velocities smaller than 50 km/s. The relative veloc¬ 
ity is small so that baryons are captured gravitationally 
by the dark matter halo and accrete to the halo. There¬ 
fore, we can roughly estimate the threshold velocity as 
the circular velocity of the dark matter halo at the initial 
redshift. 


"^cir 



Mr, 


3.85 X W Mq 


1/3 


km/s. 


( 10 ) 
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FIG. 8: The relative time difference for the collapse as a 
function of relative velocity Vbc with the halo mass fixed to 
Me ~ 4 X 10^ Mq/H and three initial density fluctuations 
Sc,i = 0.016 (black), Sc,i = 0.033 (rad) and Sc,i = 0.066 (blue) 
evaluated from turnaround time. The shaded regions show 
the standard error region from the N-body simulation. The 
solid lines are the fitting formula Eq. (HH). The dashed line 
are results from the N-body simulations with the usual pe¬ 
riodic boundary condition. The upper horizontal axis shows 
the mass converted with Vbc = 30 km/s. 


This criterion can be also obtained from the condition 
that the third term related to the relative velocity, 
dominates in the right-hand side in the second equation 
of Eq. (j2j). Therefore, when the relative velocity is larger 
than the criterion velocity, the relative motion prevents 
the collapse by the third term in Eq. ©■ Similarly, when 
the relative velocity is larger than Vdr^ the effect of the 
relative motion on the structure formation arises. The 
upper horizontal axis represents the corresponding mass 
in the case of a fixed velocity v^c = 30 km/s, which is 
converted through the v^d^c -dependence of the rel¬ 
ative motion effect. Thus one can find that the super¬ 
sonic relative motion does not affect the formation of 
dark matter halos with mass larger than Me ^ 10^ M^jh 
for Vbc = 30km/s which corresponds to the effective Jeans 
scale discussed in Sec. ttm 

In the N-body simulation results with the usual peri¬ 
odic boundary condition, the modification becomes inde¬ 
pendent of the amplitude of the relative velocity in the 
case with Vbc ^ 200 km/s shown in Fig. [HI However this 
independence is due to the artificial condition of the sim¬ 
ulations. When Vbc ^ 200 km/s, all baryons can travel a 
distance larger than the simulation box size Tbox until 
t ^ 100 and are attracted gravitationally twice by the 
collapsing dark matter sphere. Therefore, the effect of 
the relative velocity seems to be saturated in this rela¬ 
tive velocity region. On the other hand, the results of the 
N-body simulations with our boundary condition shown 
as the shaded region in Fig.|8]are not saturated. In order 


to verify the validity of the result of our N-body simula¬ 
tions, we show the collapse time and the baryon fraction 
in the limit of the large homogeneous relative velocity. 
One can easily imagine that, when the relative velocity 
is enough high {vbc = 500 km/s), which corresponds to 
very small dark matter halos (Me ^ 10^ Mq/K) in the 
case of Vbc = 30 km/s, baryons do not collapse along the 
relative velocity direction. In this limit, we can ignore 
the gravitational force from the dark matter halo on the 
baryon motion along the relative velocity direction. In 
other words, the gravitational collapsing of baryons oc¬ 
curs perpendicular to the relative velocity direction and 
does not along the parallel direction. Therefore, to eval¬ 
uate the evolution of the density fluctuations, it is useful 
to consider the motion of baryons with the relative veloc¬ 
ity in the comoving cylindrical coordinate system whose 
axis is parallel to the relative velocity motion. In this 
case, the EoM of baryons particles with the initial veloc- 
ity (V6||,f±) = (vbc,0) is given as 


d^n GSM,^ 

-w=- - w''*- 


5Mc =-7rpc[rc,i(l + 4,i) - r^] 
X max[l, [n/ref], 


( 11 ) 


where Vb is the radial component of fb from the center of 
the dark matter sphere and Vc is the radius of the over¬ 
density sphere of dark matter. Note that we ignore the 
gravitational force of the baryon fluctuation in Eq. (pT]) . 
because SMb 5Me. We solve Eqs. © and (HU numeri¬ 
cally with different initial positions fb,i chosen randomly. 
We calculate the resultant baryon density contrast in a 
collapsing spherical shell of dark matter by taking the 
average of Sb = {rbi±/Tb±)^ — 1 for baryons inside the 
shell of dark matter, because the collapse of baryons is 
cylindrical. 

Figure [9] shows the evolution of the baryon fractions 
with Me ^ 4 X IO^Mq/Zi and Vbc = 500 km/s. The red 
solid line represents the solution obtained with the peri¬ 
odic boundary condition. In order to take into account 
the periodic boundary condition, we solve Eqs. © and 
m with the assumption that the position of baryons, fb, 
is limited within the box size and baryons return into the 
box from the opposite side when they exit from one side 
of the box. For comparison, we also solve the equations 
without the boundary condition and plot the solution in 
the blue solid line. Baryons with the periodic boundary 
condition feel gravitational force stronger than without 
the boundary condition. Therefore, the collapse is faster 
with the boundary condition than without the boundary 
condition. 

Moreover, in Fig. |9l we show the results from N-body 
simulations. The red shaded region represents the stan¬ 
dard error region from the N-body simulations with the 
periodic boundary condition, while the blue shaded re¬ 
gion gives the standard error region for the N-body sim¬ 
ulations with our boundary condition which is the pe- 
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FIG. 9: The baryon fractions within the dark matter halo 
{n = 42 kpc/h) with Me ~ 4 x IO^Mq/Zi and Vhc — 500 km/s 
and solution of Eqs. © and m with or without the peri¬ 
odic boundary condition (red or blue line). The each color 
shaded region show the la dispersion of the baryon fractions 
estimated from five realizations of N-body simulations. The 
vertical dashed lines show the turnaround times from N-body 
simulations. 


FIG. 10: The collapse time with Me ~ 4 x IO^Mq/Zi and 
Vbc — 500 km/s estimated by solving Eqs. ([7j) and (ITT]) with 
boundary condition (red) and without boundary condition 
(blue). The black solid line is corresponded the solution of 
Eqs. ([7j) with (5?, = 0 at any time. Points represent estimations 
from N-body simulations with Vbe — 500 km/s. The upper 
axis presents the collapse time from Eq. o with the total 
matter density fluctuation Sm,i = Pc^c,il{pc + Pb)- 


riodic boundary condition with the position shuffling of 
the baryon particles reentering into the box. In Fig. [9l 
N-body simulations with our boundary condition is con¬ 
sistent with the numerically solution without the peri¬ 
odic boundary condition, while N-body simulations with 
the periodic boundary condition agrees with the solution 
with the periodic boundary condition. We find that, in 
the both boundary condition cases, the differences be¬ 
tween the numerical solutions and N-body simulations 
arise around i 100. This is because, as the collapse 
proceeds, the baryon density in the N-body simulations 
grows as the spherical collapse rather than the cylindrical 
one. Therefore, the baryon fraction is larger in N-body 
simulations than in the numerical calculations. 

In addition, we plot A as functions of the initial 
density fluctuations 6c,i with Me ^ 4 x 10 ^Mq/Zi and 
Vbc = 500 km/s in Fig. [TOl The red and blue solid lines 
are obtained from the numerical calculations with and 
without the periodic boundary condition, respectively. 
In the case with the periodic boundary condition, we 
overestimate the gravitational force to collapse as men¬ 
tioned above and, therefore, the delay of the collapse is 
not larger than in the case without the boundary condi¬ 
tion. For comparison, we plot the black solid line which 
represents the results with the assumption that baryons 
cannot collapse. The difference from the black solid line 
represents the contribution due the collapse of the baryon 
component. When the relative velocity is large enough, 
baryons cannot collapse to the dark matter mass shell. 
Accordingly, as the relative velocity becomes large, the 
blue solid line shifts to the black line. We also plot the re¬ 
sults of N-body simulations with the periodic boundary 
condition and our boundary condition as red and blue 


points with the standard error bars in Fig. [TOl respec¬ 
tively. As shown in Fig.O the numerical calculation with 
the periodic boundary condition agrees with N-body sim¬ 
ulations with the periodic boundary condition, while the 
numerical calculation without the boundary condition is 
consistent with N-body simulations with our boundary 
condition. We remind you that our boundary condition 
is introduced to remove the artificial distribution of the 
reentering baryon particles due to the periodic bound¬ 
ary condition in N-body simulations. We conclude that 
the saturation in A from N-body simulations with the 
periodic boundary condition is caused by this artifact. 
The results from N-body simulations with our boundary 
condition present realistic phenomena. 

Based on the results of our N-body simulations, we find 
the fitting formula of A represented as the solid lines in 
Fig. [HI The fitting formula is given by 


Ai^Mc Vbc^ 6c,i) — Asij=o{6c,i) 


{Me, Vbc, 6c,i) 
B^{Mc,Vbc,6c,i) + 1’ 


B{Mc,Vbc,6c,i) = 


Vbc 


Me 


Vnorm{6c,i) \3.Sb X 10^ Mq/H 
"^norm (^c,i) — 6 y 6 c,i, 


- 1/3 

( 12 ) 


where Asi,=o is the solution of Eq. © with 6b = 0 shown 
in Fig. [TOl and z/, ay and 6^ > 0 are the fitting param¬ 
eters. The velocity Unorm is the critical velocity for the 
collapse of baryons along the direction of the relative ve¬ 
locity. When Ubc ^ '^norm, the relative velocity is much 
larger than Ucirc even at the collapse time and baryons 
does not collapse along the direction of the relative veloc¬ 
ity as mentioned above. Since the collapse time becomes 
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long with decreasing Si and vi^c is inversely proportional 
to the scale factor, r’norm increases as Si decreases. Thus 
we conclude that the delay is controlled by two critical 
velocity Vdrc and 'r’norm- Nevertheless it can be calculated 
numerically, we use the approximated function of Asi,=o, 


Asi^=o{Sc^i) 


-1-06 

4.12 X (50;y® + 1.93’ 


(13) 


We estimate the parameters by fitting simultaneously 
Eq. ([12]) to A with using three different initial density 
fluctuations. Furthermore we perform a Fisher analysis 
and obtain the parameters as 


z/ = 2.02 ± 0.07, ay = 205 ± 16, by = 877 ± 253, (14) 



where the standard errors are estimated after marginaliz¬ 
ing over other parameters. Since we sample the data for 
three different initial conditions, the parameter by has a 
large error. However, we find that the form in Eq. m 
fits well with the N-body simulation results. 


FIG. 11: The ratio of mass function of the dark matter halo 
between with and without relative motion at four redshifts 
z = 10 (red), z = 15 (green) z = 20 (blue) and z = 30 
(magenta). 


C. Mass function 


The delay of the halo formation due to the relative ve¬ 
locity modifies the abundance of dark matter halos. In 
this section, we evaluate the modification based on the 
Press-Schechter formalism. In the Press-Schechter for¬ 
malism, the delay of the collapse is represented as the 
increase of the critical density contrast. Using the rela¬ 
tive time difference A, we can write the modified critical 
density contrast during the matter dominated era as 

^ crit c t ^c,i) — ^crit [1 T A{My^ Sy^i^] . (1^) 

The critical density contrast depends on the relative ve¬ 
locity in the region where the collapses happens. There¬ 
fore, the modified halo mass distribution can be written 
with the probability distribution function of the ampli¬ 
tude of the relative velocity at the initial time /{vijy) as 


J dvjjc f{Vbc)\j ^ 


2 p(^) ^crit('^C5 ^c,i) 


My a{My,z) 
d Inherit (Me, Vby, Sy^i) dlncr(Mc, z) 


X exp 


dMy dMy 

dcrit(Mc5 "^bci dc,i) 


2(j‘^{My,z) 


(16) 


Note that, because of the existence of the relative motion, 
the redshift of the collapse time depends not only on 
the initial density fluctuation of cold dark matter, but 
also on the halo mass My and the amplitude of relative 
velocity v^y. Therefore, in Eq. (US, the mass derivative 
term of the critical density additionally arises because the 
mass dependence of the critical density contrast affects 
the hierarchical structure formation. 


We assume that the probability distribution f{vby) fol¬ 
lows the Maxwell-Boltzmann distribution because the 
each component of relative velocity is independent and 
obey the same Gaussian distribution whose mean value 
is zero and the dispersion is cr^; 

/ 3 \ 3/2 / 3^2 \ 

f{vbc)dVbc = ^TTvl^ \ 2^) V 

where we use (7y = 28.8 km/s according to the latest 
cosmological parameters from PLANCK paper and 
we use CAMB[21| to calculate a{My^z). 

Figure [11] shows the ratio between the mass function 
with and without the relative velocity. Here we plot 
the ratio at four different redshifts, z =10, 15, 20 and 
30. These lines are evaluated by using the fitting for¬ 
mula Eqs. (jl2p for A in Eq. m- The suppression of the 
mass function due to the relative motion is more signif¬ 
icant for smaller masses region and at higher redshifts. 
At ^ = 30, although the modification A is very small 
around IO^Mq/H < My < 10^ M^jh^ the suppression of 
the mass function is not negligible. This is because such 
massive halos at high redshifts are rare objects which 
satisfy Jerit ^ o-{My) in the Press-Schechter formalism. 
Therefore, the effect of the modification A appears ex¬ 
ponentially in the Press-Schechter formalism, even if A 
is small. The ratio of the mass functions reaches the 
minimum at My ^ 10^ MqIH. On smaller scales than 
My 10^ Mq/Zz, the suppression of the mass function 
decreases because the mass derivative term of the critical 
density in Eq. m increases around My ^ 10^ Mq/H as 
shown in Fig. jS] 

Figure [TT] also shows that the mass function is sup¬ 
pressed even around the EoR {z ^ 10). This is be¬ 
cause the redshift dependence of the derivative terms 
in Eq. m is weak. Although the suppression scales 
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are consistent, the suppression of mass function around 
EoR is stronger than in other previous works. Espe¬ 
cially Ref. [13 reported that the suppression disappears 
at ^ ^ 10 in their SPH simulations. Of course, unlike 
the SPH simulations, we only simulate the gravitational 
force, ignoring the baryon physics. However, it is worth 
mentioning reasons of the difference between our calcu¬ 
lation of the mass function and previous work in term of 
the gravitational growth. The first reason for this differ¬ 
ence is that we do not include the environmental effects 
of the structure formation, e.g., an accretion of other den¬ 
sity peaks of baryons neighboring the dark matter halo. 
These effects increase the baryon fraction within a dark 
matter halo and promotes the halo formation. The sec¬ 
ond reason is that the relative motion produces halos 
derived from offsetting baryon peaks [22|. This effect 
increases simply the number of halos. It is difficult to 
include this effect in the Press-Schechter formalism. We 
will be able to estimate a number of halos originated 
from the baryon peaks by applying the our simulations 
or our semianalytical model of Eq. ©■ Einally, the ini¬ 
tial condition for matter density and velocity fields is still 
debatable in the numerical simulations with the relative 
velocity. In our simulation, we use the initial condition 
that leads to the maximum delay of the collapse time. 
Thus the halo mass function in Eig. [11] is calculated on 
the basis of the optimistic case where the baryons can 
escape most efficient from the dark matter halo. 

VI. SUMMARY 

The relative motion between baryons and dark mat¬ 
ter plays an important role, particularly, in small-scale 
structure formation at high redshifts. We have studied 
their effect on the dark mater halo formation. 

We have evaluated the delay of the dark matter halo 
collapse due to the supersonic relative motion by using 
the cosmological N-body simulation. We have found that 
the delay of the collapse becomes large for a dark matter 
halo with Me ^ 4 x 10^ Mq/Zi, when the relative veloc¬ 
ity is larger than Vdr = 57 km/s. In other words, the 
delay of the collapse happens when the relative velocity 
is larger than the typical circular velocity of the dark 
matter halo. Moreover, we have shown that the super¬ 
sonic relative motion delays the fall of baryons into the 
potential well of the dark matter halo in the context of 
the spherical collapse model. We have also evaluated the 
baryon fraction M^/M^ of the dark matter halos with 
the supersonic relative motion by the N-body simulation. 
The baryon fraction becomes smaller as the amplitude of 


the relative motion increases. We have pointed out that, 
when the relative velocity is large enough to escape from 
the potential of dark matter halos, baryons can collapse 
only along the perpendicular direction of the relative ve¬ 
locity, like the cylindrical collapse. Eurthermore we show 
the delay of the collapse time for dark matter halo by 
the relative motion depends on the initial density fluctu¬ 
ation within dark matter spheres, which determines the 
collapse time of the dark matter halo without relative 
motion. The smaller initial density fluctuation lead the 
longer time during which the supersonic relative motions 
affect the halo collapse. In consequence, the effect of 
the relative motion is more efficient on dark matter halos 
formed at later time. 

Einally we have estimated the suppression of the abun¬ 
dance of dark matter halos by supersonic relative motion 
in the context of the spherical collapse model. In the 
Press-Schechter formalism, the delay of the collapse in¬ 
creases the critical density contrast for the collapse. We 
have found the fitting formula of the critical density con¬ 
trast depending on the halo mass, the initial density fluc¬ 
tuations and, the relative velocity. Using the fitting for¬ 
mula, we have calculated the mass function of dark mat¬ 
ter halos. The relative motion decreases the mass func¬ 
tion with mass smaller than 10^ M^/h before EoR. In 
particular, the abundance of halos with Me = 10^ Mq/H 
is suppressed by 80% at z = 30 and a half at 2 ; = 10. 

The delay of the dark matter halo collapse and the de¬ 
crease of the baryon fraction in dark matter halos due 
to the relative motion can give the effect on first star 
formation and the reionization history p3l-[^. Such ef¬ 
fect could impact the cosmological signals of the EoR in- 
cluding the CMB polarization [^, the redshifted 21 cm 
lines [ 23 , [HI and these cross-correlation [H|- Moreover, 
the relative motion between dark matter and baryons in¬ 
fluences the large scale structure, e.g., baryon acoustic 
oscillation jSOl-l^ , and there is the challenging work de¬ 
tecting this effect by using the results of galaxy distribu¬ 
tion from two independent galaxy spectroscopic survey 
[H| . Based on the results of this paper, we will investi¬ 
gate the effect of the relative motion on the cosmological 
signals probed by ongoing or planned observations. 
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